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Abstract. - We investigate the behavior of an array of Bose-Einstein condensate tubes described 
by means of a Bose-Hubbard Hamiltonian. Using an anisotropic non-polynomial Schrodinger 
equation we link the macroscopic parameters in the Bose-Hubbard Hamiltonian to the ones that 
are tunable in experiments. Using a mean field approach we predict that increasing the opti- 
cal lattice strength along the direction of the tubes, the condensate can experience a reentrant 
transition between a Mott insulating phase and the superfluid one. 
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Introduction. — Since Greiner et al. first succeeded 
t storing a Bose-Einstein condensate in a two-dimensional 
array of narrow tubes using an optical lattice [1], experi- 
mental progress on low-dimensional Bose systems has been 
tremendous. To name but a few highlights, the Mott tran- 
sition was observed in three [2], one [3], and two dimen- 
sions [4,5], and in ID a Tonks- Girardeau gas has been 
realized [6,7]. Moreover, in 2D the Kosterlitz-Thoulcss 
transition was observed [8,9]. The dimensional crossover 
between three, two, and one dimensions for bosons in an 
optical lattice was studied theoretically using Tomonaga- 
Luttinger liquid (TLL) theory in Refs. [10, 11] and using 
Monte Carlo simulations in Ref. [12]. These references 
studied 2D or 3D optical lattices in which atoms can tun- 
nel easily along one Cartesian direction but not along the 
others. The system can then be described as an array of 
tubes of bosons, which may or may not be mutually phase 
coherent. As the tunneling between tubes is varied, such 
a system will undergo a transition from a 3D superfluid 
(3D SF) to a 2D Mott insulating phase (2D MI), which 
consists of a decoupled array of ID tubes. Similarly, if the 
tunneling probability along all three Cartesian directions 
are made unequal, a 2D SF state can be realized, in which 
there is superfluidity within separate 2D layers but no co- 
herence between them. These transitions are present only 
if the tubes have finite length. 



The theoretical approaches used in Refs. [10-12] are 
able to describe phase fluctuations within the tube-like 
filaments of bosons, but they do not capture any nonlin- 
ear effects due to the possible variation in the width of 
these tubes. Gross-Pitaevskii theory describes such varia- 



tions [13]; however, it does so at the expense of not being 
able to account for phase fluctuations. Keeping these lim- 
itations in mind, we offer in this Letter a description that 
is complementary to those of Refs. [10-12], using Gross- 
Pitaevskii theory in order to understand how nonlinear 
effects may affect the phase transitions in this peculiar 
type of optical lattice. Such effects become important if 
the potential barriers along the strongly coupled direction 
are so weak that each tube can be considered as a quasi- 
1D Bose-Einstein condensate and the number of bosons 
is large. It will be shown that such nonlinear effects can 
give rise to a reentrant Mott transition in the array of ID 
tubes. 

The dilute Bose gas in an external potential V ext (r) is 
described by the second quantized Hamiltonian 
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The physical setting consists of a two-dimensional optical 
lattice in the x and y plane with period d x and d y , cre- 
ating a square array of tubes which develop along the z 
direction. Moreover a weaker optical potential parallel to 
the tubes is added, so that in the present case the external 
potential is given by 
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As in Ref. [14] the many-body wavefunction of the gas 
-0(r) is rewritten as a sum of local operators 

i>(r)=J2^j v k(r)b jk , (3) 

where 4>j x j y k{*) = <i>k{x + j x d Xl y + j y d y ,z) and each 
bjk acts on the k th state of the j th tube. The (f>k are 
a complete set of wavefunctions. In the present setting 
we expect only the lowest energy state to be occupied, so 
that it is possible to drop the k index from (3). Using 
expression (3) the second quantized grand potential 



6$, ft] =H[i,,ft]-(i G N$,ft] 
can be rewritten 
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where j = (j x ,j y ), n. 
matrix elements are 



bijbj, and the inter-tube tunneling 
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+ V exk {x,y,z) <j>j a <j>j a+ i 
with a = {x, y}. The in-tube interaction energy is 

U = g j dr\^(v)\\ (7) 
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where fj.Q is a constant that fixes the number iVtot of par- 
ticles in the whole system and g = 4:7rh 2 a s /m. 

We assume that the external potential (3) along the x 
and y direction is very strong, so that the behavior of 
each tube can be described by a one-dimensional approx- 
imation of the Gross-Pitaevskii equation (cf. [13]). This 
approach is not appropriate to study the behavior of sepa- 
rate sites, since the number of particles would be too small 
to get accurate results; it is however expected to work in 
this context since each tube spans over many lattice sites. 
The non-polynomial Schrodingcr equation (NPSE) [16] 
has proven itself to be an excellent means to describe such 
many particle-systems. The main difference of this ap- 
proach compared with the TLL description [11] is that it 
is possible to consider explicitly the transverse width of 
each tube, therefore studying the nonlinear effects due to 
the finite density on the behavior of the system. The main 
assumption of this work consists in taking the groundstatc 
of the NPSE as the Wannier function <f)(x, y, z) that ap- 
pears in the expression of t, U, and fi. In this way it is 
possible to link the physical parameters that are tunable 
in a real experiment to the macroscopic ones in the Bose- 
Hubbard Hamiltonian. Notice that in general V x ^ V y , so 
an anisotropic version of the NPSE is needed. 



Anisotropic NPSE. — We suppose that the confin- 
ing potential in the transverse direction is strong enough 
so that the behavior of each tube is well described consid- 
ering a harmonic external potential 



Vtubc(r) 
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Therefore the wavefunction of each tube is factorized as 



(r) 



where 
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We note that the approximation of the Wannier functions 
by Gaussians is known to introduce quite large quanti- 
tative errors when computing the tunneling matrix ele- 
ments [13]. Nevertheless, since the Gross-Pitaevskii equa- 
tion does not allow an easy treatment using Bloch theory, 
it is very hard to improve on the Gaussian approximation. 

Following the steps in Ref. [16] an equation for the longi- 
tudinal part of each tube is derived. The energy functional 
for the Bose-Einstcin condensate in each tube is 
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Substituting (12) into (14) it is possible to carry out the 
integration along x and y. Dropping the terms propor- 
tional to rj'(z) [15] and minimizing with respect to r\ and 
/* the following equations are obtained: 
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Vy (z) = a y (l + 7 |/(z)| 2 ) 3 , 



where a y = \fhj (rnu> y ) and 
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Notice that Ref. (15) reduces to the previously obtained 
NPSE [16] in the case A = 1. 
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Bose-Hubbard Parameters. Using (12) as an 
ansatz for the condensate wavefunction it is possible to 
express t x , t y , fj,, and U in terms of integrals of f(z) and 
rj(z). The expression for t y is 
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Fig. 1: Left panel: behavior of t/U as a function of V z for 
Wccii = 1,2,3,4,5, while V x — V y = 26.4_Z5r; for a system 
of length L the parameter is scaled according to (24). Right 
panel: MF order parameter (6) against V z for the case n cc n = 4 
and L = 6 (2.55/im). 
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Consistently with the NPSE approximation the terms 
proportional to r]'(z) and r]'(z) 2 have been omitted. The 
expression for t x is obtained from t y switching the indices 
x and y, and changing A to A -1 . The integral for U gives 
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Switching the indices x and y leaves the value of U invari- 
ant, so that the expression for the self- interaction energy 
of each tube is consistent with the symmetry of the prob- 
lem. 

In order to make explicit calculations we consider a sys- 
tem of 87 Rb atoms, whose atomic interaction is repulsive 
(a s = 5.77 nm). The lattice is a simple cubic one, with 
d z = d y = d z = 425nm. Equation (15) is solved for f(z) 
and r](z) with periodic boundary conditions over a single 



cell by means of a self-consistent approach. The length of 
the tube affects the numerical results through the normal- 
ization of the function f(z), so that using periodic bound- 
ary conditions over L cells the value of t a (a = {x, y}) is 
left invariant and U scales down as 1/L, i.e. ti, = t while 



(23) 



where t l and Ul denote the values of t and U for a system 
of tubes with length L. Notice that 



ta \ ta 



(24) 



Denoting by n ce ii the number of particles in each longitu- 
dinal well, then the number of particles n tu be in each tube 
is 

^tubo = L n coll . (25) 

At first we focus on an isotropic case, i.e. a situation 
in which V x = V y = 26.4 En, where £r is the recoil en- 
ergy for a wavelength of 780 nm. Defining t = t x = t y , 
in Fig. 1 (left panel) there appears the behavior of t/U 
while varying V z from to 10 -Er, for n ce u = 1,2,3,4,5. 
According to Fig. 1, for n ce ii = 4 the ratio t/U decreases 
until V z « 7 Er, then it begins to rise. It is known [17, 18] 
that the Mott insulating phase appears as a set of "lobes" 
in the plane fi/U-t/U, each one corresponding to a pre- 
cise number of particles n tu be for each tube; outside of 
the 2D MI zone the system moves along lines of constant 
^tube- The initial decrease and subsequent decrease of 
t/U therefore means that for an appropriate length of the 
tubes the system can cross the 3D SF - 2D MI transition 
twice. In order to give an estimate of the critical poten- 
tial strength a mean field (MF) approach is employed (see 
for example [18,19]). The MF order parameter is given 
by the expectation value of the destruction operator (b). 
Figure 1 (right panel) plots the behavior of this parameter 
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Fig. 2: Behavior of i (left panel) and U (right panel) as a 
function of V z for n oe li = 4, while V x = V y = 26.4_Er. 



for n ce ii = 4 and L = 6 (2.55/im). The bounded interval in 
which the superfluid phase vanishes is a clear consequence 
of the nonlinear behavior of the system. 

In order to understand the physical meaning of this phe- 
nomenon, the plot in Fig. 2 shows the behavior of t and U 
against V z separately for the case n cc ii = 4. Increasing the 
longitudinal optical lattice f(z) becomes narrower, thus 
raising the value of U (22), while at the same time the 
wavefunction widens in the radial direction, increasing the 
tunneling rate t (19). At first U rises faster than t but for 
t ks 7 Er this relation is reversed. The observed nonlinear 
effect is therefore the result of a competition between t 
and U. 

So far it has been considered the isotropic case, in 
which t x is identical to ty. That is no longer true in an 
anisotropic setting, i.e., one in which V x ^ V y . Figure 3 
(left panel) shows the behavior of t x and t y upon vary- 
ing V y , with an occupation number n co n = 4, keeping V x 
and V z fixed at 26.4 Er and 5 Er respectively. It is shown 
that while t x varies very little for different values of V y , 
t y decreases exponentially. In the right panel of Fig. 3 is 
plotted in detail t x against V y for the same physical pa- 
rameters, showing that t x has a minimum for V y ~ 35 Er. 
These results suggest that upon increasing V y the boson 
gas enters a 2D SF phase, in which the system is orga- 
nized in superfluid layers along the y direction. Moreover 
the behavior of t x suggests that - if all the parameters are 
carefully tuned - these layers could experience a reentrant 
transition between the 2D SF phase and the 2D MI one, 
in which each tube is isolated from the others. 

The asymmetric case provides another interesting phe- 
nomenon. Figure 4 shows the situation in which V x = 
26.4 E R , and rc ce u = 4, for V y = 26.45 E R (left panel) and 
V y = 26.5 Er (right panel). Increasing V y , the curves of t x 
and t y become more and more separated; in particular the 
curve for t y shifts downwards, and its minimum moves to 
the left. This behavior suggests a curious effect in which - 
upon increasing V z - the system may move from a 3D SF 
phase to a 2D SF one, from this to a 2D MI phase and 
then back again to the 2D SF. 
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Fig. 3: Left panel: plot of t x /U and t y /U as a function of V y , 
while V x = 26.4 .Er, V z = 5 -Br, and n ce ii = 4. Right panel: the 
behavior of t x /U is plotted in more detail for the same physical 
parameters. 



Conclusions. — In this Letter we have studied the be- 
havior of a 2D array of strongly elongated tubes of bosons 
in an optical lattice, employing the Bose-Hubbard Hamil- 
tonian in order to describe the 3D SF - 2D MI transi- 
tion. Using the NPSE (15) to compute the groundstate in 
each tube it was possible to link the parameter t and U 
of the Bose-Hubbard Hamiltonian to the physical param- 
eters of the atom gas. Using this approach an observable 
nonlinear effect is found: we predict that in an array of 
87 Rb tubes 2.55 /Ltm long, with d x = d y = d z = 425 nm, 
V x = V y = 26.4 Er, and an occupation number of 4 atoms 
for each cell, the system goes through the 3D SF - 2D MI 
transition twice while increasing V z . The system itself is 
found to be in the insulating phase for V z between ~ 5 Er 
and ~ 8£fl. In addition, in the anisotropic case, our re- 
sults suggest that the nonlinear behavior of the system 
should cause a reentrant transition between the 2D SF 
phase, in which the gas is organized in superfluid layers, 
and the 2D MI phase, where each tube acts independently. 

A real experiment will be complicated by the fact that 
not all the tubes would be equally long and the occupation 
number fluctuates. Moreover, due to the approximations 
we have made, the most severe being the assumption of a 
Gaussian profile for the basis functions and the neglect of 
quantum fluctuations in the phase, the exact numbers are 
expected to differ from our predictions. We argue that 
the qualitative analysis should hold, i.e., an experiment 
should in certain parameter regimes show a dip in the 
phase coherence while varying V z , since this feature de- 
pends crucially only on the fact that the tunneling t and 
on-site interaction energy U exhibit different functional 
dependencies on the width of the tube-like condensates. 
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Fig. 4: Behavior of t x /U and t y /U as a function of V z while 
V x = 26.4 E R and n ce ii = 4; left panel: V y = 26.45 £ R , right 
panel: V H = 26.5 En- 
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